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Aims. To understand the fundamental physical processes important for the evolution of solar rotation and 
Q . distribution of chemical species, we provide theoretical predictions for particle mixing and momentum 
CZJ . transport in the stably stratified tachocline. 

ci ■ 

• • Methods. By envisioning that turbulence is driven in the tachocline, we compute the amplitude of turbulent 
flow, turbulent particle diffusivities, and eddy viscosity, by incorporating the effect of a strong radial 

X 

» I differential rotation and stable stratification. We identify the different roles that the shear flow and stable 



stratification play in turbulence regulation and transport. 

Results. Particle transport is found to be severely quenched due to stable stratification as well as radial 
differential rotation, especially in the radial direction with an effectively more efficient horizontal transport. 
The eddy viscosity is shown to become negative for parameter values typical of the tachocline, suggesting 
that turbulence in the stably stratified tachocline leads to a non-uniform radial differential rotation. Similar 
results also hold in the radiative interiors of stars, in general. 
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1. Introduction 

Since the formation of the radiative core, which marked the beginning of its journey on the 
main sequence, the sun has slowed down significantly due to the loss of angular momentum 
from its surface (e.g. see Stix 1989; Schatzman 1993). The angular momentum transport must 
have been very efficient during its spin-down in order for the sun to have rotational profile as 
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observed today (see, e.g. Charbonneau et al 1998). Vigorous turbulence in the convection zone 
and possibly thermal wind (Miesch, Brun & Toomre 2006) can readily provide a mechanism for 
efficient radial momentum transport, thereby eradicating radial differential rotation therein. Such 
turbulent transport is however considered to be absent in the stably stratified radiative interior, 
which has also spun-down during the solar evolution, presently rotating uniformly at the rate 
roughly the same as the mean average rotation rate on the solar surface. Whichever mechanism 
is responsible for momentum transport in the interior (which itself is an important problem), 
it should be closely related to the transport in the tachocline through which the surface spin- 
down is communicated to the interior. Transport in the tachocline also plays a crucial role in 
the overall mixing of light elements (lithium, beryllium, etc) (see e.g. Schatzman 1993; Brun, 
Turck-Chieze, & Zahn 1999), thereby determining the level of their surface abundances on the 
sun. Therefore, it is essential to understand physical mechanisms for transport in the tachocline 
and then to formulate a consistent theory starting from first principles based on those processes. 
This is particularly true since virtually all the previous theoretical modelling heavily relies on a 
simple parameterization of transport process, which is then adjusted to obtain the agreement with 
observations. 

We have initiated the development of a consistent theory of turbulent transport in the 
tachocline in the previous papers, by taking into account the crucial effect of the large-scale shear 
flows, provided by a strong radial differential rotation (Kim 2005) as well as latitudinal differen- 
tial rotation (Leprovost & Kim 2006). By envisioning that the tachocline is perturbed externally 
[e.g. by plumes penetrating from the convection zone above (e.g. see Gilman 2000; Brummell, 
Clune & Toomre 2002; Rogers & Glazmaier 2005)], or by instabilities (e.g. see Watson 1981; 
Charbonneau, Dikpati & Gilman 1999), we demonstrated how turbulence level and transport 
are reduced via shearing in a non-trivial manner (see also Burrell 1997; Kim & Diamond 2003; 
Kim 2004;2006) in a simplified three dimensional (3D) hydrodynamic turbulence. In particular, 
turbulent transport of chemical species and angular momentum are shown to become strongly 
anisotropic with effectively much more efficient transport in the horizontal (latitudinal) direction 
than the vertical (radial) direction due to shear stabilization by strong radial shear. The resulting 
anisotropic momentum transport was shown to reinforce a strong radial shear (Leprovost & Kim 
2006), with a positive feedback on the confinement of the tachocline (Spiegel & Zahn 1992) 
while chemical species are predicted to have latitudinal dependent mixing due to the variation 
of radial shear (Kim 2005). Furthermore, the results indicate that the turbulence regulation by a 
shear flow (i.e. differential rotation) leads to weak turbulence and mixing in turbulent tachocline. 

The purpose of this paper is to investigate how a stable stratification in the tachocline modi- 
fies the predictions obtained in these studies. We again envision that the tachocline is turbulent, 
driven by a forcing as in Kim (2005). In the presence of a stable stratification, the turbulence 
in the tachocline is no longer completely random as the stable stratification provides a restoring 
force against radial displacement of fluid elements, supporting the propagation of internal gravity 
waves (Lighthill 1978). These waves tend to increase the memory of otherwise random turbulent 
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fluid motion and can reduce the overall transport due to turbulence. We shall show that the tur- 
bulent transport due to shear stabilization found in the previous studies is further enhanced in the 
presence of stable stratification (gravity waves). Thus, gravity waves acting together with shear 
stabilization can lead to a weak mixing in the tachocline, as required for the surface depletion of 
lithiums (e.g. Pinsonneault et al. 1989). 

It is important to contrast our approach to those adopted in most previous works, which 
focus on the momentum transport by gravity waves themselves through dissipative processes 
(e.g. Plumb 1978; Kim and MacGregor 2001;2003; Talon, Kumar & Zahn 2002). For instance, 
gravity waves are considered to be generated in the convection zone (Press 1981) and deposit 
their momentum into background shear flow via radiative damping as they propagate through 
the tachocline and the interior (Kim and MacGregor 2001;2003). Relying crucially on radiative 
damping, the momentum transport by gravity waves would occur on a long time scale. [Recall 
that the transport by waves requires molecular dissipation, and is thus a slow process.] Similarly, 
weak mixing due to damped gravity waves was also suggested as a mechanism for a modestly 
enhanced mixing of light elements (lithium) (Press & Rybicki 1981; Garcia Lopez & Spruit 
1991). In these works, a certain level of background turbulence could not be ruled out and was 
invoked to provide an enhanced viscosity for the evolution of a mean flow (Kim and MacGregor 
2001;2003) and the gravity waves (Charbonnel & Talon 2005; Talon, Kumar & Zahn 2002). The 
value of the effective viscosity is often arbitrarily taken to be a positive constant, being much 
larger than molecular viscosity, not affected by shear flow nor by gravity waves, or it is simply 
parameterized. 

In this paper, we shall treat turbulence and gravity waves on an equal footing and consis- 
tently compute the values of turbulent eddy viscosity and particle diffusivity, by incorporating 
the effect of a shear flow (provided by a radial differential rotation), and identify different roles 
that gravity waves and shear flow play in turbulent transport. Specifically, we shall show that 
unlike shear flows, which reduce both turbulence level and transport, a stable stratification can 
suppresses turbulent transport without much effect on turbulence level. We shall further demon- 
strate that the stratification favors a negative eddy viscosity and thus tends to sharpen the radial 
gradient of large-scale shear flows rather than smoothing it out. This tendency was also found in 
Kim & MacGregor (2001;2003). We note that the elucidation of the effects of gravity waves is 
essential for understanding the momentum transport in the radiative interior as they have often 
been advocated as a mechanism to explain a uniform rotation in that region. Stratified turbu- 
lence with a shear flow is also important for the transport in radiative interiors and/or envelops 
of stars, in general, as well as in geophysical systems. In particular, it has actively been studied 
in geophysical systems (e.g. see Jacobitz, Sarkar & van Atta 1997; Stacey, Monismith, & Burau 
1999), where stable stratification was shown to inhibit turbulent transport in the direction of a 
background density gradient, leading to the two dimensional (2D) turbulence. 

The remainder of the paper is structured as follows. We first investigate the effect of gravity 
waves on turbulence in Sect. 2. In Sect. 3, we incorporate the effect of strong radial differen- 
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tial rotation and study how the gravity waves modify the overall turbulent transport. Section 4 
contains the conclusion and discussions. 



2. Internal Gravity Waves 



U (x) 




Fig. 1. The configuration of our model. 

To simplify the problem, we shall consider incompressible fluid with local cartesian coor- 
dinates x,y, and z for radial, azimuthal, and latitudinal directions, respectively (see Fig. 1), and 
use Boussinesq approximation to capture the effect of stratification. We assume that the fluid is 
stirred by a forcing on small scales, giving rise to fluctuations. In the absence of stratification, 
this fluid forcing will drive turbulence on small scales and maintain it at the level at which the in- 
jected energy is balanced by dissipation in the system. As the stratification increases, the forcing 
will generate not only random turbulent motion but also coherent (gravity) waves. Alternatively, 
some of turbulent motion will be turned into packets of gravity waves. Therefore, we consider 
both turbulence and gravity waves as fluctuations on scales much smaller than those associ- 
ated with mean density or mean background shear flows. Specifically, we express the total mass 
density p = po(x) + p' where po(x) and p' are the background and fluctuating mass densities, 
respectively; the total velocity u = Uo + v where Uo and v are a large-scale shear flow due to 
differential rotation and small-scale fluctuations; the total particle density of chemical elements 
n = «o(x) + n' where «o(x) and n' are mean and fluctuating components. Then, the main govern- 
ing equations for fluctuations v, p', and n', involving both turbulence and gravity waves, are as 
follows (see, e.g., Kim & MacGregor 2003; Moffatt 1978): 



(d, + Uo ■ V)v = -Vp - gp'x + vV 2 v + f , 

V-v = 0, 

(<9 f + U -V)p' = P^ Vx +pV 2 p', 
g 

(dt+Uo-vV = -d x n v x + DV 2 n' . 



(1) 
(2) 
(3) 
(4) 



Here, v, p and D are molecular viscosity, thermal diffusivity, and particle diffusivity, respectively; 
f in Eq. (]]) is the small-scale forcing driving turbulence; p = po(x = 0) is the constant background 
density [measured at the bottom of the convection zone (e.g. see Kim and MacGregor 2003)], 
and N = ^-g(d x po + pglc\)fp is the Brunt-Vaisala frequency, where c s is the sound speed. Note 
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that the typical values of v, D, p, and N in the tachocline are 10 2 cm 2 s _1 ,10 2 cm 2 s _1 , 10 7 cm 2 s _1 , 
and 3 x 10~ 3 s _1 , respectively. 

In this section, we first ignore a background shear flow and study how the turbulence and 
transport are modified due to stable stratification in the tachocline. To obtain the overall turbu- 
lent transport, we solve coupled equations dTJ-([4j> in terms of Fourier transform for fluctuating 
quantities (/>': 

(f>'(x, t) = — J— r f cPk </>(k, t) exp {i(k x x + k v y + k z z)} . (5) 
(2ny J 

Equations (Q]i-® then give us an equation for p = e vk ~'p in the form 

d„p + (p. - v)k 2 ~p + N 2 p = —e vkh h i . (6) 

g 

Here, N 2 = yN 2 /(y+a 2 ). Since N » pk 2 is valid on a broad range of reasonable scales / » 10 4 ~ 
10 5 cm in the tachocline, we can find the solutions to leading order in (vk 2 /N <SC) pk 2 /N <K 1 as: 



V x ~ - 

y 



-!— 5- f dt'GJt,l!)cosN(t-t'jkx{t'), (7) 
+ a 2 J 

--Vx~- f dt'G(t,t')h 2 (t'), (8) 
7 7 J 

-—v x +- f dt'G(t,t')h 2 {f) , (9) 
7 7 J 

N 



' f , 
dt 

+ a 2 )J 



G„(f,f')sin/V(r-f')MO- (10) 



Here, G(tj') = e ->* 2 ('-'') and G^tj') = e -(M^(f-f)/2 the Green's functions; a = kjk y , 
P = k z /k v , and y — 1 + 1 ; h\ and in Eqs. © and (TTb-ifTOb are the forcing terms which are 
related to fi as 

h = (l+p 2 )f x -af y -apf z , (11) 
h 2 = -pf y + f z . (12) 

Equations (IT}-© show that the vertical motion v x , involving G^ only, is always subject to radia- 
tive damping p while the horizontal motions v v and v x , containing both G M and G, are much less 
affected by p.. 

For simplicity, we assume the forcing to be homogeneous in space with a short correlation 
time Tf. 

(fi(ki,h)fj(k2, ? 2 )> = T f (2n) 3 6(t, - ? 2 )<J(k, + k 2 )^ i7 (k 2 ) . (13) 

Note that in the case where the forcing f due to plumes induces a stronger turbulence towards 
the bottom of the convection zone, f becomes inhomogeneous with the power spectrum ifr^ = 
t//ij(x, k) depending on the radial coordinate x. 

A long but straightforward algebra by using Eqs. (T7b — (fl~3b then gives us the following results 
on turbulence level and particle diffusivities defined by (n'v,) = -D'jdjUQ. 
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0n(k) 



) 3 2(y + a 2 ) 2 (ft + v)k 



<V?) ~ Tf 



a 0n (k) 



J (2tt; 
r d 3 i 

J (2n 

r d 3 i 
T f r dH 

N 2 J (27) 3 47(7 + a 2 ) 

r d- 
Tf J (^ 

r d 
Tf J (^ 



1 



TP 

TP 
3/ 
7T 
3/ 



1 



y6 2 022(k) 1 



27 2 (7 + a 2 ) 2 (// + v)k 2 
a 2 /? 2 0n(k) 1 
27 2 (7 + a 2 ) 2 (jj. + v)k 2 

0n(k) 



2wfc 2 



"22 



(k) 1 



2wfc 2 



)V 

3 /v 



a 2 0n(k) | /3 2 22 (k) 
4^(7 + « 2 ) 2 2v(v + Z)) 2 ^ 
fl 2 ^ 2 n (k) 22 (k) 



(14) 
(15) 
(16) 
(17) 

(18) 
(19) 



AN\y + a 2 ) 2 2v(v + W 
Here, we assumed spatial symmetries (pij(k y ) = (pij(-k y ) and <pij(k y ) = 0y(— k inside the 
integrals in Eqs. dT4b — (fT9b is the wavenumber of the forcing; 0y(k) (i = 1,2) is the power 
spectrum of the forcing defined by: 



(hi(kuti)hj(k2, f 2 )> = Tfdnfsih - «a)tf(ki + k 2 )0, 7 (k 2 ) . 



(20) 



In this paper, we shall consider an incompressible forcing (V ■ f = 0) for simplicity, in which case 
<f>ij in Eq. ( 1201) is related to 0y in Eq. dT3l > as: 

k 4 1 
^11 = ^4 011 ,012 = Tj(k 2 k x k z 4r u +k 2 k 2 H ip n ), 



022 = Tj(klk 2 z i// U + ^033 + 2k x k z k 2 H iff 13 ) . 



(21) 



Note that 0n and 22 can signify strong radial forcing (e.g. by plumes) and horizontal forcing 
[e.g. due to the instability of the latitudinal shear in the tachocline (see also Kim 2005)]. For 
instance, in the case of strong radial forcing by plumes (0n » 033) with ku » k x (see Fig. 2 in 
Leprovost & Kim 2006), 0n will dominate over 22 . Note further that in the 2D limit with k, = 
and 033 = 0, 22 vanishes. 

We discuss some of the important implications of Eqs. (TT4b — <fT9l). First, all three components 
of fluctuating velocity amplitude in Eqs. (TT4b-(fT6b are independent of N, clearly showing that the 
amplitude of the turbulent flow is not influenced by stable stratification in the case of the forcing 
with a short correlation time. The exact level of fluctuations is determined by the characteristics 
of the external forcing (0,y). We examine this in the simple 2D limit where k, = (7 = 1) and 
033 = 0. In this limit, the substitution of 22 = in Eqs. (fT4li-(fT6]i gives <v 2 )/<v 2 ) ~ 1/a 2 = k 2 lk 2 x 
and (v 2 ) = 0. This is an expected result for the 2D incompressible fluid in x — y domain. If the 
forcing f contains only gravity modes, the wave number of the forcing would satisfy the local 
dispersion relation k x = +k v -sjN 2 /(aj - Uok y ) 2 - 1 (here, lj is the frequency of the gravity waves). 
Furthermore, if these gravity modes are generated by the overturning of fluids in the convection 
zone (e.g. Press 1981) with Uq ~ 0, they would have a strong power at low frequencies a> <SC 
N, thereby giving k x /k y ~ N/a> » 1, and thus (v 2 )/(v 2 ) <k 1. Note that this is the situation 
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normally considered in the previous works (e.g. Kim & MacGregor 2001;2003; Talon, Kumar 
& Zahn 2002), where the main focus was on the momentum deposition by such gravity waves 
due to radiative damping after entering the tachocline from the bottom of the convection zone. In 
contrast, in this paper, the forcing f is not restricted to gravity modes, but is taken to be general, 
including any form of perturbation with arbitrary values of k x /k y . For instance, in the case of a 
strong radial forcing due to plumes with 0n » 022 and a = k x /k y <k 1 (see Fig 2 in Leprovost 
& Kim 2006), (v 2 )/(v 2 > » 1 with a stronger radial fluctuation than horizontal one. Furthermore, 
Eqs. (fT4T> -<fl~6b show that the level of anisotropy measured by the ratio of the turbulence amplitude 
depends on v and fi. For instance, for an isotropic forcing with <f>n ~ <p22, (v 2 .)/^ 2 ) °c v/jj. «: 1 
becomes small as the radiative damping yiz increases. This is because a large fi tends to decrease 
(vertical) turbulence level by introducing large (thermal) dissipation. Thus, without a shear flow, a 
stronger horizontal turbulence (level) than vertical one can be caused by large thermal diffusivity 
(but not by stratification) in the case of a temporally short correlation forcing. Note that the 
anisotropy in turbulence level could be related to the Peclet number Pe = vl/fi as (v 2 )/(v?) Pe. 
Here, v and I are the characteristic velocity and length scale of the forcing (which is fixed). 

Second, turbulent transport in Eqs. (ITTb — (TT~9b are strongly affected by stratification as N 2 
increases, in contrast to fluctuation levels in Eqs. (fT4b— (fT6b. Since N » fik 2 on reasonable scales 
/ » 10 4 ~ 10 5 cm, the vertical (radial) mixing is severely quenched as the stratification increases 
(i.e. as N~ 2 ), in qualitative agreement with Brun, Turck-chieze & Zahn (1999). In comparison, 
the part of the horizontal (latitudinal) mixing due to the radial forcing tf> n is reduced proportional 
to N~ 2 while the one due to the horizontal forcing 022 is independent of N. The comparison of 
these two contributions (or alternatively D xx and D yy ) in the case of <f>u ~ 022 gives us a cut- 
off scale l c ~ y/v/N above which vertical (radial) mixing is strongly reduced compared to the 
horizontal (latitudinal) mixing. That is, in the presence of both radial and horizontal forcings 
of comparable strength, stable stratification mainly reduces the vertical transport without much 
effect on horizontal transport on scales / > l c . For parameter values typical of the tachocline 
v ~ 10 cirrs~ and N ~ 3 x 10 -3 s -1 , l c ~ 10 cm. Thus, the stratification is very likely to play 
an important role over a broad range of physically reasonable scales. In the limit of strong radial 
forcing with <p\\ » 022 and a ~ 0, the stratification influences both the radial and horizontal 
transports to the same degree while the incompressibility renders Dj'/Dj ~ Dj'/Df ~ 1/a 2 ~ 
(v 2 )/(v 2 ) » 1, with an effectively more efficient radial transport. It is interesting to compare 
our result D xx cc N~ 2 with the mixing due to radiatively damped waves (e.g. Garcia Lopez 
and Spruit 1991; Talon et al 2002). For instance, Garcia Lopez and Spruit (1991) estimated the 
vertical mixing due to gravity waves to be proportional to n/N 2 . While the reduction in vertical 
mixing for large is in qualitative agreement with our results (oc Af~ 2 ), the increase in vertical 
mixing in Garcia Lopez and Spruit (1991) is due to the fact that damped waves are necessary 
for wave transport. Finally, we note that without a shear flow, momentum transport vanishes (i.e. 
(v x v y ) = (v x v z ) = 0) for an isotropic forcing. Scaling of turbulence amplitude, turbulent viscosity 
(vr), and turbulent diffusivity Dj are summarised in Table 1 for an isotropic forcing. 
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A large thermal diffusivity (ji = 10 5 v) in the tachocline is often considered to reduce the 
stabilizing effect of stable stratification via the weakening of the buoyancy restoring force. To 
highlight this effect, it is illuminating to consider the extreme limit of strong thermal diffusion 
where density fluctuation becomes stationary with d,p' = N 2 vjg - fik 2 p' = in Eq. (f3]). In this 
limit, by following a similar analysis as previously (with D = v), we can easily obtain the vertical 
and horizontal particle diffusivities as follows: 



where a = jN 2 /fik 2 (y + a 2 ). Equation (l22l shows that the vertical mixing D xx oc fi 2 /N 4 oz Pe~ 2 
decreases for large N while increasing for large yu. This is because the reduction in the vertical 
mixing due to buoyancy force is weaken by a strong radiative damping id. The comparison of Eqs. 
(|22"1 i-(|2"3|> further shows that the reduction in vertical mixing relative to horizontal mixing is given 
by a factor of (vk 2 ) 2 /a 2 for an isotropic forcing and is thus weaker than that in the case of weak 
radiative damping fik 2 <k N [see Eqs. (fT7]i-(fT9]ll. Furthermore, this reduction appears on scales 
larger than the critical scale l Cil = (jj.v) 1 ^ 4 /N 1 ^ 2 ~ 10l c . Here, l c = y/v/N is the critical scale in the 
case of fik 2 «: N. These results thus show that a strong thermal diffusion weakens the buoyancy 
effect and makes the effect of stratification become important on larger scales, compared to the 
case of a weak thermal diffusion. This result is thus consistent with the expectation employed in 
previous works. 

To summarize, a stable stratification can dramatically quench turbulent transport with a more 
effective mixing in the horizontal directions orthogonal to the background density gradient. It 
does not however affect the amplitude of the turbulent flow, the ratio of which is found to depend 
only on v/fi (for a temporally short-correlated forcing). 

3. Consistent theory 

The results in Sec. 2 showed that for a temporally short correlated forcing, a stable stratifica- 
tion reduces turbulent transport only, leading to anisotropic turbulent transport, without much 
effect on turbulence level. In this section, we study how these results are modified by a stable 
background shear flow (differential rotation in the tachocline). In particular, we will show that a 
shear flow not only inhibits the vertical mixing further, enhancing the anisotropic transport, but 
also reduces turbulence levels anisotropically, thereby leading to effectively stronger horizonal 
turbulence. For simplicity, we ignore the latitudinal differential rotation compared with the radial 
differential rotation since it is weaker in the tachocline due to thin tachocline (h < 0.03 ~ 0.05 of 
the solar radius R). The inclusion of the latitudinal differential rotation would introduce a small 
correction term in our results. For instance, for turbulence amplitude, this correction term is of 
order (h/R) 2 <s 1, as shown in Leprovost and Kim (2006). Note that the latitudinal shear is 
crucial for non-vanishing horizontal momentum transport in Leprovost & Kim (2006). Again, 




(22) 



(23) 
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,-2/3 
bv 





Table 1. Scaling of turbulence amplitude, turbulent viscosity (vr), and turbulent diffusivity Dj 
for an isotropic forcing with <p\\ ~ <p22 in the case N » fj.k 2 and D ~ v «: yu. ^ = fik 2 ,/Jl <sc 1, 
£ v = vk 2 JJ[ <K 1, and £o = Dk 2 JJ[ <K 1 are small parameters representing strong shear limit. 
The second (M = and N + 0), the third (M * and N = 0), and the fourth (M * and N + 0) 
columns contain the results for stratified unsheared case in Sec. 2, for unstratified sheared case 
(Kim 2005), and for stratified sheared case in Sec. 3, respectively. 



we envision that turbulence is maintained in the tachocline by an external forcing while gravity 
waves are excited due to this external forcing in the stably stratified tachocline. We shall then 
compute the overall turbulent transport consistently by taking into account the interaction among 
turbulence, shear flow and gravity waves, instead of simply assuming a (large) constant value of 
turbulent viscosity for mean shear flow (and gravity waves). Note that this treatment is essential 
when there is no clear scale separation between gravity waves and turbulence, in which case 
turbulence cannot be considered to give an enhanced value of viscosity for gravity waves (c.f. 
Charbonnel & Talon 2005). 

For the evolution of fluctuations, we approximate the radial differential rotation by a linear 
shear flow with Uo = -xtty to keep the analysis tractable. Here, is the shearing rate which 
we assume to be positive without loss of generality. As done in previous papers (Kim 2005; 
Leprovost & Kim 2006), to capture the effect of shearing due to radial differential rotation {J[ ~ 
3 x 10~ 5 s _1 for the tachocline) non-perturbatively, we use the special Fourier transform for 
fluctuating quantities <p'\ 

<t>'(x, t) = ——; \ d 3 kip(k, t) exp {i(k x (t)x + k y y + Lz)} . (24) 
(2ny J 

Here, k x = k x (t) is the time dependent [unlike constant k x in Eq. ©], satisfying an eikonal 
equation 



d,k x (t) = k y M . 



(25) 
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Equation d25l l implies that k x linearly increases in time as k x (f) = k x (0) + k y J[t, manifesting the 
main effect of shearing by a shear flow Uo(x)y, i.e., generation of fine scales in the x direction 
due to tilting and distortion of fluid eddies (e.g. see Burrell 1997; Kim 2004; Kim 2005). The 
efficient generation of fine scales by shearing leads to the breakup of eddies and enhancement of 
the overall dissipation, thereby reducing turbulence amplitude and transport (e.g. see Kim 2005). 
A similar effect by a shear flow is expected to persist for a more realistic radial shear [e.g. for an 
error function used in helioseismic inversions (Kosovichev 1996; Corbard et al 1999)] since the 
basic mechanism of shearing (e.g. see Fig. 2 in Kim 2005) is the same regardless of the details 
of the profile of radial shear. The efficiency of the shearing could depend on the details of the 
profile, possibly leading to a slightly different scaling. 

It is interesting to note that a gravity wave with an initially positive value of k x (0) can change 
its sign after the time interval k x (Q)/k v ^l (for k y > 0) due to shearing. Since the local (ra- 
dial) group velocity of gravity waves is given by v gx = -k 2 k x N 2 /k 4 (u> - Uok y ) (e.g., see Kim 
& MacGregor 2003), the gravity wave thus alters its propagation direction as k x flips its sign. 
Therefore, a gravity wave which initially propagates downward to the interior from the convec- 
tion zone with a negative vertical group velocity (i.e. v gx < 0) can propagate upwards when the 
vertical group velocity becomes positive (y gx > 0) due to shearing. 

For parameter values typical of the tachocline N ~ 3 x 10~ 3 s and J[ ~ 3 x 10~ 6 s , 
Ri = N 2 /3l 2 » Ri c =i 1/4, satisfying the stability criterion (Lighthill 1978). We thus assume 
that the radial shear flow is stable with large value of Ri = N 2 /^R 2 » 1 in the remainder of the 
paper (see also Schatzman, Zahn, & Morel 2000). Note that the buildup of chemical composition 
gradient (the so-called fj. gradient in the solar interior [e.g. see Michaud and Zahn (1998) and 
references therein] would further increase the values of N and Ri, making the radial shear flow 
more stable, although this effect could be counteracted by a radiative damping, as shown in Sec 
2. The shearing rate J[ ~ 3 x 10~ 6 s _1 due to this radial differential rotation is larger than the 
dissipation rate due to radiative damping /uk 2 on a broad range of scales l(= l/k y ) > 10 6 ~ 10 7 
cm ~ 10~ 3 Hq, where Ho is the pressure scale height ~ 6 x 10 9 cm. Thus, we focus on the strong 
shear limit in the following by using ^ = fj.k 2 /Jl « 1 as a small parameter. 

For Ri » 1 and v ~ D <K fi, a long but straightforward algebra can give us the solutions 
to Eqs. (Q~|)-(|4|), as shown in Appendix A. By using these solutions and the correlation functions 
of forcing defined in Eq. (1201 . we obtain the following results for turbulence level and transport 
coefficient defined by (n'v,) = -Djdjno for i = 1, 2 and 3, and momentum flux (v x v y ) = —vt8 x Uq 
in the strong shear limit ^ = [ik 2 jy{ <K 1 (see Appendix A for details): 




(27) 



(26) 



(28) 
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°-*Jw^ &G " (29) 



77 r A 

r 2 J (2tt)V 



''ii(k)r(|) p / 3 V | yg^ 22 (k)r(f ) / 3 \" 



2N 2 y^y + a 2 V \&/ ^ 2 \ 2 & 



3 Jt 0n(k) 

v r - 



) 3 y(y + a 2 ) 



1 

+ 



J?l 2 12Af 2 (y + a 2 ) 



(30) 
(31) 



Here, f v = v£ 2 A?l « 1, £ D = Dfe 2 /^ « 1, £ M = ntf/M « 1, G = |r(l/3) (3/2^ v ) I/3 

1 / \2/3 

and Gi = (2/3) ( . Note again that ^ « 1 is valid on a broad range of scales l(= 
1/ky) > 10 6 ~ 10 7 cm in the tachocline and also that ^ « 1 guarantees that ^73 ~ £ v «: 1 since 
D ~ v <& /j. The spectrum 0y in Eqs. d26l>— <f3TT> are given by Eqs. (120) , which are related to 
power spectrum of forcing tfr/j in Eq. dT3~b in the incompressible case. 

Equations (f2~6b-(f3~Tl> reveal the following interesting features. Turbulence levels given in Eqs. 
d26t -(l28b are again independent of stratification, similarly to the case without a shear flow [Eqs. 
(TT~4b — <TT~6bl while they are reduced for strong shear 31. This indicates that waves and shear flows 
play different roles in turbulence regulation - waves do not necessarily quench fluctuation levels 
while shear flows can reduce them through enhanced dissipation via shearing. The turbulence 
regulation by shearing gives us horizontal velocity fluctuations in Eqs. d27b and d28l i which are 
effectively higher than vertical one in Eq. d2"o*i i. While a similar tendency was also found in 
the absence of gravity waves (Kim 2005), the exact value of the ratio of vertical to horizontal 
turbulence levels is not the same. For example, (v 2 )/(v 2 ) oc (| ln^,|) _1 for a strong radial forcing 
0ii with 022 = while (v 2 )/(v 2 ) oc ^ 3 for an isotropic forcing with 0n ~ 022- Note that 
%v <K ^ <K 1 are small parameters in our problem, representing the strong shear limit. These 
results are to be compared with (v 2 )/(v 2 ) oc in the unstratified medium (Kim 2005). The 
results for the isotropic forcing in various cases are summarised and compared in Table 1 . 

Transport properties in Eqs. d29ll-(l3"TT> however exhibit a very different behaviour, with 
both vertical and horizontal mixing being inhibited in a non-trivial manner by strong stratifi- 
cation as well as by shearing. First, the vertical transport is reduced as D X j oc ^d^^N~ 2 oc 
D[T 2 ^3{~ l ^N~ 2 , becoming small as either stratification or shearing increases. Note that the de- 
crease in Dj for large N agrees with Miesch (2003). Interestingly, Dj oc jjT 2 ^ oc Pe 2 ^ decreases 
as the radiative damping fi increases. This is because large radiative damping increases thermal 
dissipation, thereby mainly inhibiting the vertical mixing, as noted in Sec. 2. Thus, compared to 
the case with N = where D" oc 3\T 2 (Kim 2005), the vertical mixing is much more quenched 
by a factor of % D ^ 2/3 (3l/N) 2 « f D ^ 2/3 = (D/fi)^J 3 since 31/N < 1 and D/// ~ 10~ 5 

in the tachocline (see also Table 1). This clearly shows that shear flow (orthogonal to radial 
density gradient), stable stratification, and radiative damping can all inhibit the radial transport. 

Second, Z>^ is less affected by stratification since it involves the two parts - the one from 0n 
is proportional to ^o^^N^ 2 while the other from 022 is proportional to yr 2 £J 2 ^ 3 , independent 
of N. In the simplest case of a strong radial forcing with 022 = 0, D^/D^ ~ £/ 3 <k 1, indepen- 
dent of N. This is similar to the case without stratification where D"/Dj ~ t^J 3 (Kim 2005). In 
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the general case where 0n + and 022 + 0, exhibits a non-trivial, complex interplay among 
stratification and shear flow in determining the overall transport. To appreciate this, we compare 
Eq. < f30l > with the result obtained without shear flow ( fT8l ) to find that shear flow enhances the 
contribution from 0n by a factor ^d^ 4 ^ = (Df/u)^ 1 ^ while it reduces the part from <f>22 by a 
factor of (v/,?[) 2 £; 2/3 cc ^ /3 <k 1. Since D/fi ~ 10~ 5 <s 1 and «: 1, the horizontal transport 
driven by a strong radial forcing 0n can be either inhibited or enhanced by a shear flow (due to 
0i i) depending on the parameter values. On the other hand, the horizontal transport due to 022 
is reduced by the shear flow by a factor proportional to <K 1 . The end result can easily be 
shown to be the enhancement of the effect of stratification. To see this, we compare the two con- 
tributions from 0ii and 022 to Dj and obtain a characteristic scale 1* = (D/ fi) 2 (v/N) ' 2 (tt/N) ^ 2 
above which the stratification is important, being mainly responsible for quenching transport. 
Since I, is smaller than l c = (v/N) 1 ^ 2 obtained without a shear flow, the effect of stratification 
becomes important compared to the case without shear flows. It is also worth comparing the con- 
tributions to Dj from 0n and 022 separately to those in the case of N = in Kim (2005), where 
Dj cc J{r 2 ^y 2li from both 0n and 022- That is, the contribution from 0n is further reduced by a 
factor of (yi/A^) 2 ^y 3 (v//i) 4 ^ 3 by stratification, while the one from 022 is not affected. 

We now examine how the stratification affects the anisotropy in transport. Equations ( f29T > and 
(f30b for a strong radial forcing with 0n » 022 and a ~ gives D^/Dj ~ ^ 3 «: 1, independent 
of N. This is analogous to the result D^/D^ ~ if 2 / 3 <K 1 obtained in the unstratified case (Kim 
2005). In this case, the anisotropy in the transport is solely caused by shear stabilization. For an 
isotropic forcing with 0n ~ 022 and a ~ 1, Dj/D^ ~ (JM/ N) 2 {v I n) 2 ^ % v «: 1. In other words, 
the anisotropy in turbulent transport depends on stratification, shear, and v/fi, as clearly shown 
in Table 1. Since (Jl/N) 2 «: 1 and v//i « 1 in the tachocline, Dj X /D^ is much smaller than 
£ 2 ^ 3 obtained in the unstratified case (Kim 2005). That is, the anisotropy in transport is further 
enhanced due to stratification. We emphasize that the anisotropy in turbulent transport is much 
stronger than that in turbulence amplitude, discussed previously (see also Table 1). This result 
also shows the reduction in Z>™ for large fj. and again highlights the importance of the radiative 
damping in reducing the vertical transport, thereby increasing the anisotropy in the transport. 

Finally, Eq. OTb demonstrates one of the most important effects of a stable stratification, 
which is to drive a system away from a uniform rotation with a negative eddy viscosity. Recall 
that in the absence of stratification the eddy viscosity is positive in 3D while negative in 2D limit 
(e.g. see Kim 2005). In contrast, the eddy viscosity in Eq. OTT ) is negative in both 2D (J3 — and 
y — 1) and 3D cases. This behavior was also found in recent numerical simulation by Miesch 
(2003) of a stably stratified turbulence with an imposed shear driven by penetrative convection, 
who found anti-diffusive radial and diffusive latitudinal momentum transport, thereby offering 
a mechanism for a proper transition from latitudinal differential rotation in the convection zone 
to solid body rotation in the radiative interior. Anti-diffusive momentum transport is a generic 
feature of a strongly stratified medium (i.e. a geophysical system). It is interesting to note that a 
negative viscosity was also found in the previous work on momentum transport due to radiatively 
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damped gravity waves (e.g. Kim and MacGregor 2001;2003). In addition, the result QTl ) shows 
that the (anti-diffusive) momentum transport becomes less efficient for strong stratification (large 
AO, in agreement with Miesch (2003). 

To summarise, our results show that the shearing effect by radial differential rotation together 
with gravity waves is an important mechanism for turbulence regulation, leading to a weak tur- 
bulent transport and anisotropic turbulence and transport in the tachocline. Furthermore, we have 
consistently derived the values of turbulence level, particle mixing and momentum transport 
starting from the first principle, clearly identifying the different roles of gravity waves and shear- 
ing in transport. For instance, in comparison with Chaboyer and Zahn (1992) [or Spiegel and 
Zahn (1992)] which start with the assumption of a strong horizontal mixing, we identified the 
source of such an anisotropic turbulence. In particular, we have made a clear distinction between 
the anisotropy in the turbulence level and turbulent transport. 

4. Discussion and conclusions 

We have studied turbulent transport in the stably stratified tachocline with a strong radial dif- 
ferential rotation, when turbulence is driven and maintained by a forcing (e.g. due to plumes 
penetrating from the convection zone or due to instability). We have assumed that both turbu- 
lence and gravity waves are on small scales (with no clear scale separation between the two) 
and treated the interaction among gravity waves, turbulence and shear flow consistently. Unlike 
a shear flow which regulates both turbulence level and transport, a stable stratification is shown 
to mainly inhibit turbulent transport, leading to a further reduction in transport compared to 
the unstratified case (Kim 2005). Specifically, for parameter values typical of the tachocline 
(N/^l » 1), particle transport due to a strong radial forcing <p u (with 022 = 0) is reduced as 
^o^ 2 ^ N~ 2 , becoming much smaller than horizontal transport (oc ^ D ^ A ^N^ 2 ) by a factor of 
<sc 1. Here, £ v = vk 2 /Jl <K I, g D = Dk 2 /Jl <sc 1, and ^ = fik^/M <K & are the small 
parameters characterizing a strong shear limit (see the main text for more details). Note that in 
this case D^/Dj oc ^ 3 <g 1 depends only on shearing but not on stratification, indicating that 
the anisotropy in particle transport is mainly governed by radial differential rotation. A similar 
scaling (oc % 2 ^ « 1) was also found in the unstratified case (Kim 2005). However, in the case of 
an isotropic forcing with <p\\ ~ <p22 and k x /k y ~ 1, the horizontal mixing is much less reduced 
with Dj ~ 3{T 2 tfy 2 ^ (with no effect of stratification), leading to a stronger anisotropy in trans- 
port with D^/D^ ~ (Jl/N) 2 (v/iJ.) 2/3 ^y <k £; /3 <K 1. Note that the anisotropy becomes stronger 
for larger fi. 

Furthermore, the vertical momentum transport was shown to be anti-diffusive with a nega- 
tive eddy viscosity. That is, small-scale turbulence influenced by gravity waves accentuates the 
gradient in a radial differential rotation rather than makes it uniform. This is similar to the ten- 
dency obtained in the case of momentum deposition by gravity waves due to radiative damping 
(Kim & MacGregor 2001;2003). The sharpening of the gradient of the radial shear due to the 
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negative viscosity could eventually lead to time variation in the tachocline (similarly to Kim & 
MacGregor 2001) or instability (e.g., see Petrovay 2003), causing a rapid radial mixing and thus 
reducing the anisotropy. 

Even in the stably stratified radiative interiors of stars as well as in the tachocline, background 
turbulence has often been assumed to be present to enhance the value of effective viscosity. While 
the clarification of the source of turbulence responsible for such an enhanced eddy viscosity 
is an interesting problem, our results demonstrate how this residual turbulence interacts with 
gravity waves and shear flow, providing the prediction for the values of eddy viscosity as well 
as the particle diffusivity which depend on physical quantities like 3i and N. In particular, the 
radial particle transport D™ ~ ^ D ^ 2/3 N^ 2 indicates that turbulent transport of particles can be 
inhibited due to stable stratification, shear flow and large radiative damping. This finding can 
have interesting implications for the surface depletion of lithium in the Sun and other stars (i.e., 
Pinsonneault 1997), and will be studied in a future publication. Furthermore, if a similar physical 
process operates in the bulk of the radiative interior of the sun, a negative viscosity that we obtain 
implies that a radial differential rotation which is created during its spin-down would not be 
eliminated. Note however that Charbonnel & Talon (2005) have shown an efficient momentum 
transport in the solar radiative interior due to the cumulative effect of large scale meridional flow, 
shear instability and gravity waves. 

The tachocline is believed to possess a strong toroidal magnetic fields of strength 10 4 ~ 10 5 
G, which could have an important influence on turbulent transport in that region. This is an impor- 
tant problem since the presence of a weak poloidal magnetic field in the radiative interior together 
with a strong toroidal magnetic field in the tachocline (acting as a boundary layer between the 
radiative interior and convection zone) could offer a mechanism for a uniform rotation in the in- 
terior as well as for the tachocline confinement (e.g. Riidiger & Kichatinov 1996; MacGregor & 
Charbonneauu 1997; Gough & Mclntyre 1998) [see however Brun & Zahn (2006) for a negative 
result on this scenario]. In this case, the values of effective diffusivity of magnetic fields as well 
as eddy viscosity play a crucial role in determining the thickness of the tachocline. Therefore, 
a consistent computation of magnetic diffusivity and eddy viscosity with their dependences on 
physical quantities such as the strength of magnetic fields, the Brunt-Vaisala frequency, and 
shearing rate, would be of primary interest (Kim and Leprovost 2007). It is an interesting ques- 
tion, in general, whether a negative eddy viscosity, favored in a stably stratified medium, remains 
as a robust feature in the presence of magnetic field in view of the forward energy cascade (i.e. 
positive eddy viscosity) in MHD turbulence. Furthermore, composition gradients (discussed in 
Sec. 3) and meridional flows would contribute to the transport in the tachocline. In particular, 
meridional flows are expected to enhance the radial transport of chemical species by advection 
as discussed in Kim (2005) although this enhancement would be reduced for stronger horizon- 
tal turbulent [e.g., see Kim (2005) and Chaboyer and Zahh (1992)]. These issues are however 
outside the scope of this paper and will be addressed in future publication. 
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Appendix A: 



In this Appendix, we show how to derive Eqs. (f26b— (f3Tb. To this end, we use Uo(x)y = —xMy in 
Eqs. ([T)-©, and introduce the transforms w and w for w = V;, n, p, hj and / as follows: 



w = vvexp[v(^/3fc^?t + fc|*)L 
w = wexp[p(k 3 x /3k y Jl + k 2 H t)/2]. 



(A.l) 



where k 2 H - k 2 + k\. We then use a new time variable r = + J[t to rewrite Eqs. (HJ-© as: 



Md T v x = -irkyp - gp + f x , 

&d T Vy - MV X = -ikyp + fy , 

J\d T v z = -ik z p + f z , 

= TV X + V y + — V z , 



m[d T + f»(y + T 2 )] P = — v,, 



N 



&[d T + t D (y + T 2 )]fi = -(d x n )v 



(A.2) 
(A.3) 
(A.4) 
(A.5) 

(A.6) 
(A.7) 



For a strong shear limit where p = pk 2 /J[ <K 1, Eqs. ( IA.2b — ( lA76l can be combined to yield the 
equation for p as 



BAiy + AdM + ^-p 



N 



ft 2 r g rt 



■Ai(t), 



(A.8) 



to leading order in ^ = [xkZ/M <K 1. Here, h\{T) = yh x — rh y - t/3Ii z ; y — 1 +/3 2 and fi = k z /k y . In 
the limit of N 2 /$l 2 » 1, Eq. ( IA.8b can be solved with the following solutions valid up to 0(N 2 ): 

.72 



P(T) = — 2-777 2T7 -Sin[^(T)-¥>(Ti)]/Zi(Ti), 

^(t + t 2 ) 5 / 4 J C 7(Ti)(r+T?) 1 



+N-Jy + t 2 vj(t) cos [^(t) - ^(ir)] 



--sin [<p(r) - <p(Ti)] 



P 2 *, , P 

p(r) 

ytr yft 



Vy(T) = P*(t) + — p(r) — (iTlG(T,Tl)fe(Tl) 

y Jvy 
TP* P 



v z (t) 



V,(t) - 75-p(T) + — £/TlG(T,Tl)fe(Tl). 

y yiy yyt 



(A.9) 



Here, again /3 = kjk y , y = 1 + p 2 ; N 2 = N 2 /gJi, N = N 2 y/Jl 2 - G(t,ti) = 
exp{(^/2 - £,)[(yr + r 3 /3) - (yr, + rf/3)]}; f v = v^/^l; 

/ll = (1 + y6 2 )/l x - n%y - Tlp% z , 
ft 2 = -fitly + , 



, s Na , Jy + T 2 + t 1 r 
p(T) = — In- 



2 Jy + T 2 - t 8iV Jy + r 2 ' 
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w(t) = a ■ 



J 2_ 



a = 1 - 



(A. 10) 



8N 



By following a similar algebra used in Kim (2005) and by using the correlation function of the 
forcing given in Eq. (l20i >. we can obtain the following correlation fictions (v,V/) to leading orders 
in 1/N 2 «c 1 and^ = jj.k 2 /Jl <sc 1: 



<v 2 > 



2 N t / r <^ 3 ^ t 
<Vy) ~ # J w 

2 TV T (PkdT 



-f„Q(T, fl ) I'll 



(k) 



2Vr + f7r 2 (r + T 2 ) 3/2 



+ „-2f„e(r, fl ) ^ 2 022(k) 



7 



e -f,6(T. a ) frlQO ^ , J -2f„Q(r, fl ) 022(k) 

y 2 



Wy) — 7TTT 



2Vr + a I ^ 2 ^ + r2)3/2 

3 /Wt 0n(k) e -^ Q(T ' a h 



Vr+^^(fl) 2 (r + T 2 ) 5 / 2 

^ sin [<p(t) - <p(a)] + N -\jy + T 2 vj(f) cos [fir) - <p(a)] 



(A.ll) 



Here, a = k x /k y , Q{j,a) = (yr + t 3 /3) - (ya + a 3 /3); special symmetries iffij(ky) = ij/iji—ky) and 
tyijikz) = >frij(-k z ) were used. We then compute r integrals in Eq. ( IA. lib in the strong shear limit 
& <sc ^ «=; 1 and use <v,v v > = -v^ x d x U = vfM to obtain Eqs. (|26]l-(|28]l and d3~Tl 
Next, to compute particle transport, we integrate Eq. ( IA.7b in time to obtain 



n(k(f), = -3,n I dt^k^g (k, f; kj, ti)e- DG(Vl) V;(ki, x, fj) . 



(A. 12) 



Here, g(f, fi) = £ dt'[k 2 (f) + k 2 H ] = [k\ - k] x ]/3k y Jl + k\(t - t x ); k\ = k 2 y + k\ is the amplitude 
of wave number in the horizontal plane; k 2 - k 2 H + k 2 ; g is the Green's function given by 

g(k, t; kuh) = 6(k y - k ly )6(k z - k lz )S [k x - k u - k ly (t - ti)&] . (A.13) 

A similar analysis using Eqs. dA.9l ), dA. 13b , (fTT])-(fT3l), and d20"l ) and (nv,-> = -Djjdjn gives 



/ 

D yy = f d i kdT lE 

T (2tt) 3 M 2 J y 2 



2(2tt) 3 J1 2 N 2 



dhdr^ ( ^l Q( ™\y + ^ 2 ., 



i(k)r 2 ^j + t 2 



B -1f v Q(T,a)n2 



2N 2 y[y~ 



'/3^22(t - a)\y + t 2 ) 



+ a L 



(A. 14) 



Here, %r, = Dk 2 /J[ ~ £ v <k 1. Finally, the evaluation of r integrals in Eq. ( I A. 141 ) gives us Eqs. 
(|29]l-(|30]l in the main text. 
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